require(PBSadmb) # tool for reading admb output (R_report in code)
setwd("d:/home/_mymods/msvpa/msm2011/")
msm1<-readList("orig_diet_r.rep")
msm2<-readList("tindex_diet_r.rep")
msm3<-readList("tbiomass_diet_r.rep")
source("R/Plot_Sel.R")
source("R/Plot_M.R")
source("R/agefits.R")
names(msm3)
Plot_Sel(msm1)
Plot_M(msm1)
Plot_M(msm2)
Plot_M(msm3)

par(mfcol=c(1,1))
names(msm1)
plot(1979:2510,msm1$Biomass_1,ylim=c(0,4e7))
plot(1979:2510,msm1$Biomass_2,ylim=c(0,3e6),ylab="Pacific cod",xlab="Year",type="l",lwd=3,col="red")
plot(msm3$Biomass_2)
plot(msm1$Biomass_3)

names(msm1)
Plot_M(msm1,spp=2)
msm3$pollock_consumed

plot(1979:2010,msm3$M2_2[,1]+msm1$M1_2[1],ylab="M1 + M2",xlab="Year",ylim=c(0,1))
lines(1979:2010,msm3$M2_2[,1]+msm2$M1_2[1],col="red")
lines(1979:2010,msm3$M2_2[,1]+msm1$M1_2[1],col="black")
lines(1979:2010,msm3$M2_2[,2]+msm1$M1_2[2],col="red")
lines(1979:2010,msm3$M2_2[,3]+msm1$M1_2[3],col="blue")
lines(1979:2010,msm3$M2_2[,5]+msm1$M1_2[5],col="green")
lines(1979:2010,msm3$M2_2[,7]+msm1$M1_2[7],col="yellow")


msm1$M2_1[,4]+msm1$M1_1[4]
names(msm1)

Nage3_atf=read.table("nage3_atf.dat")
Nage3_pcod=read.table("nage3_pcod.dat")
Nage3_poll=read.table("nage3_poll.dat")

msm1$srv_age_like_2



msm1$obs_catch_hat_1
plot(1979:2011, msm1$N_1[,3],ylim=c(0,1.9e7),type="l",ylab="Pollock N age 3",xlab="year")
points(1979:2010, 1000*Nage3_poll,pch=19)
legend(1998,1.9e7,c("2010 Assessment","MSM"),lty=c(1,-1),pch=c(-1,19))
  
plot(1979:2011, msm1$N_2[,3],ylim=c(0,1.e6),type="l",ylab="Pcod N age 3",xlab="year")
points(1979:2010, Nage3_Pcod,pch=19)
legend(1995,1e6,c("2010 Assessment","MSM"),lty=c(1,-1),pch=c(-1,19))

plot(1979:2011, msm1$N_3[,3],ylim=c(0,.6e6),type="l",ylab="Arrowtooth flounder N age 3",xlab="year")
points(1979:2010, Nage3_atf,pch=19)
legend(1995,.6e6,c("2010 Assessment","MSM"),lty=c(1,-1),pch=c(-1,19))
 
Plot_Bhat(msm3)
AgeFits(msm3)
AgeFitsSrv(msm3)
AgeFits(msm3,sppname="Pacific cod",spp=2)
AgeFits(msm3,sppname="Arrowtooth",spp=3)

AgeFitsSrv(sppname="Pacific cod",spp=2)
AgeFitsSrv(sppname="Arrowtooth",spp=3)

detach(dat)

# package that does error bars easily...
require(plotrix)

source("~/_mymods/crab/crab-model/R/tools.r")
msmfit=read.fit("msm") # Separate utility to read in par, cor, std files etc

  source("/home/rhome/common.r")
  tab1=cbind(cm1$N)
  colnames(tab1)<-c('New crabs', 'Pre-recruits', 'Recruits')
  rownames(tab1)<-cm1$yrs_pf
  xtab(tab1,caption="Table 1.  Numbers at stage estimated from the model.",dec=c(0,1,1,1),cornername='Year', file="d:/t.html",width=.5)
  tab2=cbind(1:6,cm1$L)
  rownames(tab2)<-c("Avg. wt in catch","Trawl survey", "Pot survey","Trawl survey composition","Pot survey composition","Fishery composition")
  xtab(tab2,caption="Table 2.  Likelihood components.",dec=c(0,3),cornername='Component', file="d:/t2.html",width=.6)
  for (i in 1:5)
  {
    cat(runs[,i],file="cm.ctl")
    system("cm -nox")
    nam<-paste("cm",i,sep="")
    tmp=readList("cm_r.rep")
    assign(nam,tmp)
  }
for (i in 1:5)
{
  nam<-paste("cm",i,sep="")
  print(nam)
  cm.figs(dat=nam)
}
cm.figs(cm1)
cm.figs(cm2)
cm.figs(cm3)
cm.figs(cm4)
cm.figs(cm5)

for (i in 1:5) print(get(paste("cm",i,sep=""))$M98)
detach(dat)
rm(dat)

dev.off()

  cm1$yrs_pf
  tab3=cbind(cm1$L,cm2$L,cm3$L,cm4$L,cm5$L)
  colnames(tab3)=c("Model 1","Model 2","Model 3","Model 4","Model 5")
  rownames(tab3)<-c("Avg. wt in catch","Trawl survey", "Pot survey","Trawl survey composition","Pot survey composition","Fishery composition")
  xtab(tab3,caption="Table 3.  Likelihood components.",dec=c(3,3,3,3,3),cornername='Component', file="d:/t3.html",width=.6)
?matrix
tab4=matrix(nrow=35,ncol=5)
tab4[1:34,]=cbind(cm1$recs,cm2$recs,cm3$recs,cm4$recs,cm5$recs)
dim(tab4)
tab4[35,]=0
tab4[35,]=c(mean(cm1$recs),mean(cm2$recs),mean(cm3$recs),mean(cm4$recs),mean(cm5$recs))
colnames(tab4)=c("Model 1","Model 2","Model 3","Model 4","Model 5")
rownames(tab4)<-c(cm1$yrs_pf,"Mean")
rownames(tab4[35:35,])<-"Mean"
rownames(tab4)
xtab(tab4,caption="Table 4.  Recruitment estimates.",dec=c(0,0,0,0,0),cornername='Year', file="d:/t4.html",width=.6)
tab4

c("Avg. wt in catch","Trawl survey", "Pot survey","Trawl survey composition","Pot survey composition","Fishery composition")

dim(cm1$N)

  

